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GETTING  THE  GOOD  BOUNCE:  TECHNIQUES  FOR  EFFICIENT 
MONTE  CARLO  ANALYSIS  OF  COMPLEX  REACTING  FLOWS 


ABSTRACT 


The  direct  simulation  Monte  Carlo  method  as  it  is 
applied  to  complex  reacting  flows  is  discussed  at 
length.  Relations  are  presented  for  the  efficient 
treatment  of  energy  dependent  cross  sections,  multiple 
species  and  chemical  reactions  in  multi-dimensional 
flows.  Special  emphasis  is  placed  on  recent  advance¬ 
ments  in  the  method  developed  at  Spectral  Sciences, 
including  the  use  of  dynamically  adjustable  weighting 
factors  and  global  time  counters. 


DTIC  QUALITY  INSPECTED  3 


Icoosslon  For  | 

:  STIS  GRA&I 
.  DTIC  TAB 

J  Unannounced 

I  Justification. 

0^ 

□ 

□ 

i _ _ _ _ _ _ _ 

1 

1  Rv 

|_pi St ribut 1 on/ 

i  AvallabilltJ 

Codes 

jAvall  and/or 
Dlat  I  Special 


□  □ 


1.  OVERVIEW  OF  THE  DIRECT  SIMULATION  MONTE  CARLO  METHOD 

The  direct  simulation  Monte  Carlo  method,  as  pioneered  by 
G.  A.  Bird, provides  a  powerful  technique  for  the  simulation 
of  real  gas  flows.  It  bridges  the  gap  between  continuum  and 
free  molecular  flow,  retaining  validity  in  either  extreme.  It 
can  be  used  to  describe  complex  mixtures,  including  effects  of 
chemical  reactions,  heat  conduction,  viscosity  and  diffusion 
for  flows  in  three  dimensions.  To  date,  it  is  the  only  approach 
which  can  demonstrate  these  abilities  for  general  flow 
configurations. 

The  basic  calculational  technique  is  well  described  by  its 
originator  in  Ref.  1,  to  which  frequent  reference  will  be  made. 
The  present  purpose  is  to  describe  how  the  technique  is  imple¬ 
mented  at  Spectral  Sciences,  Inc.  (SSI),  with  special  emphasis 
on  extensions  developed  at  SSI  and  elsewhere  after  the  publica¬ 
tion  of  Ref.  1.  Elementary  concepts  and  relations  which  are 
essential  to  a  coherent  explanation  are  included  here  for 
clarity. 

The  direct  simulation  Monte  Carlo  method  involves  storing  a 

discrete  number  of  molecules  (via  their  velocities,  positions 

and  other  pertinent  information)  in  a  computer.  The  solution 

region  is  broken  up  into  a  number  of  separate  cells,  and  the 

solution  is  stepped  forward  in  time  in  a  two  stage  process. 

First,  the  molecules  are  advanced  along  their  trajectories  by 

an  amount  appropriate  to  their  velocity  and  a  time  increment 

At  .  In  this  first  stage  some  molecules  will  leave  the  solution 
m 

region  and  some  will  be  introduced  as  determined  by  the  boundary 
conditions  for  a  particular  problem.  The  second  stage  is  to 
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simulate  collisions  in  each  cell  appropriate  to  At^  so  that 
collision  frequencies  are  properly  simulated.  A  basic  hypothesi 
of  the  method  is  that  if  the  time  step  is  made  small  enough  the 
processes  of  translations  and  collisions  can  be  uncoupled  in  thi 
manner. 

Periodically,  the  solution  is  sampled  by  accumulating  sta¬ 
tistical  sums  of  niimber  densities,  velocities  and  other  basic 
properties.  The  solution  is  run  repeatedly  until  statistical 
deviations  are  reduced  to  a  desired  limit,  and  then  physically 
meaningful  output  quantities  are  computed  from  the  statistical 
sums.  The  number  of  molecules  represented  is  typically  a  few 
thousand  at  a  time,  which  is  vastly  fewer  than  the  number  of 
molecules  occurring  in  virtually  all  real  flows.  Hence,  the 
construction  of  a  dynamically  similar  flow  to  be  simulated  in 
the  computer  is  an  essential  feature  of  the  method. 


2.  GAS  MODEL  AND  EQUILIBRIUM  PROPERTIES 


2.1  Preliminary  Equilibrium  Gas  Relations 

For  most  problems  of  interest  there  is  a  far  field  equil¬ 
ibrium  state  whose  properties  are  of  relevance  to  the  problem  to 
be  solved.  Frequently  length  and  velocity  scales  for  the  prob¬ 
lem  are  obtained  from  the  far  field  state  and  used  to  nondimen- 
sionalize  the  internal  code  variables.  Even  when  the  far  field 
state  is  not  used  for  scaling  purposes,  it  still  provides  an 
important  comparison  case. 

For  a  rest  gas  in  equilibrium  the  normalized  distribution 

function  for  the  relative  speed,  c  ,  between  molecules  of  spe- 

.  (2) 

cies  i  and  molecules  of  species  j  is  given  by' 


where 


(1) 


a 


ij 


JO 


2R_T 

0  <» 


(2) 


and  u .  . 
ij 


is  the  reduced  mass  of  the  pair; 


i.e. 


m^m 


j 


^i 


(3) 


with  m^  and  m^  representing  the  masses  of  the  two  species.  In 
these  relations,  is  the  far  field  temperature  and  is  the 
universal  gas  constant.  (Rq  is  used  instead  of  Boltzmann's 
constant  since  the  molecular  masses  will  be  consistently 
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represented  in  atomic  mass  units  rather  than  grams.)  The  avail¬ 
able  translational  collision  energy  between  two  molecules,  E^, 
is  given  by 


E  =  i  y . .c^ 
c  2  r 


(4) 


2 . 2  Ana''ytical  Form  of  the  Collision  Cross  Section 

Whenever  the  direct  simulation  Monte  Carlo  method  is  applied, 
it  is  necessary  to  make  tradeoffs  between  accuracy  and  simplicity 
in  molecular  models.  It  does  no  good  to  use  a  complex  molecular 
potential  surface  and  then  find  that  reasonable  computer  run 
times  result  in  very  large  statistical  fluctuations  in  the  out¬ 
put.  Since  the  final  output  will  reflect  errors  in  the  statis¬ 
tics  as  well  as  errors  in  the  models,  there  is  a  strong  impetus 
to  use  models  which  contain  the  essential  physics,  but  which  can 

be  applied  in  a  computationally  efficient  manner.  The  current 

(3) 

state-of-the-art  is  the  Variable-Hard-Sphere  (VHS)  model. 

In  this  model  molecules  have  a  collision  cross  section  which 

varies  as  an  inverse  power  of  the  available  collision  energy. 

Hence,  if  a. .  is  the  collision  cross  section  for  collisions  of 
ID 

species  i  with  species  j ,  then  is  given  by  a  relation  of  the 

form 


=  A. .  e”^ 
ID  c 


(5) 


where  A^^  is  a  constant  coefficient.  It  follows  that  the  effec¬ 
tive  diameter  for  molecules  of  species  i,  d^,  is  implicitly 
defined  as  a  function  of  available  collision  energy  by  the 
relation 


a 


ii 


A.  .E 


IX 


-0) 

c 


(6) 
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can  be  determined  from  a  reference  cross  section  and  velocity 

via 

=  [“ii  •  <’> 

If  a  reference  cross  section  is  given  for  a  reference  temperature 
rather  than  a  reference  velocity/  then  the  usual  choice  for  the 
reference  velocity  is  that  velocity  which  has  a  collision  energy 
equal  to  the  mean  collision  energy  occuring  in  collisions  at  the 
reference  temperature.  Mathematically,  this  is  equivalent  to 


c 


c 


r  11 


a.  . 

r  11 


(8) 


where  the  bars  over  the  quantities  indicate  averages  taken  over 
the  distribution  function  given  in  Eq,  (1)  evaluated  for  m^  =  m^ 
and  Equation  (8)  can  be  simplified  to  give 


4(2  - 


m. 

1 


(9) 


For  simulations  involving  a  large  number  of  species,  refer¬ 
ence  cross  sections  are  frequently  not  available  for  all  possible 
collision  pairs.  In  this  case  it  is  possible  to  specify  A^^^  for 
self  collisions  only,  and  then  use  Eq.  (6)  to  get  a  molecular 
diameter  as  a  function  of  collision  energy.  Then,  applying  the 
relation 

a..  =  1T  [(d.  +  d  .)/2]  (10) 

the  coefficient  in  Eq.  (5)  for  interspecie  collisions  is  given  by 


6 


(11) 


For  the  internal  workings  of  a  Monte  Carlo  code,  it  is 
usually  rore  convenient  to  express  the  collision  cross  section 
as  a  function  of  the  relative  collision  velocity  rather  than  the 
collision  energy.  This  is  simply  achieved  via  the  relation 


a 


ij 


where 


-2oj 

r 


(12) 


(13) 


The  parameter  o)  can  be  related  to  n,  the  exponent  of  distance 

.  (3 ) 

in  an  inverse  power  intermolecular  force  law  via  the  relation 


=  2/(n  -  1) 


(14) 


Hence,  hard  sphere  molecules  (for  which  n  goes  to  infinity)  are 

represented  by  oj  equal  to  zero.  There  is  a  substantial  body  of 

evidence,  however,  that  the  effective  size  of  molecules  does 

indeed  decrease  with  increasing  collision  energy,  so  a  positive 

value  of  u)  is  usually  a  better  choice.  w  can  be  determined  from 

molecular  beam  data,  or  from  its  macroscopic  implications.  For 

example,  if  s  is  the  temperature  exponent  for  the  coefficient  of 

(3) 

viscosity,  then  it  can  be  shown  that 

s  =  w  +  1/2  ,  (15) 

so  a  measurement  of  the  temperature  dependence  of  the  viscosity 
coefficient  serves  as  an  indirect  determination  of  w. 


In  order  to  incorporate  the  model  for  internal  energy  trans¬ 
fer  to  be  discussed  in  Section  4,  it  is  necessary  that  cu  be 
assumed  the  same  for  all  interactions.  This  represents  one  of 
the  major  restrictions  in  the  current  state  of  modeling. 
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Although  the  sizes  of  molecules  are  allowed  to  vary  in  the 
VHS  molecular  model  in  deciding  whether  or  not  a  collision  is  to 
occur,  when  a  collision  does  occur  the  post  collision  velocity 
components  are  computed  as  if  it  were  a  hard  sphere  collision 
(see  Section  4).  This  results  in  a  substantial  computational 
simplification  and  yet  retains  good  agreement  with  the  macroscopic 
predictions  of  the  more  exact  model. (See  Ref.  1  for  a  dis¬ 
cussion  of  molecular  scattering  for  general  power  law  potentials.) 


2 . 3  Equilibrium  Reference  Properties  for  a  Multi-Component  Gas 

One  advantage  of  the  VHS  model  is  that  the  molecules  have  a 
well  defined  cross  section,  so  it  is  possible  to  define  a  mean 
free  path  without  putting  limitations  on  the  minimum  deflection 
angle  that  is  considered.  As  is  the  general  case  for  multi- 
component  gases,  however,  each  component  has  its  own  mean  free 
path,  and  the  overall  mean  free  path  for  the  mixture  must  be 
defined  as  a  weighted  average  of  the  mean  free  paths  of  the  indi¬ 
vidual  species.  The  somewhat  cumbersome  relations  required  to 
calculate  the  overall  mean  free  path  are  given  here.  It  should 
be  noted  that  the  mean  free  path  is  calculated  only  once  for  a 
given  problem,  so  the  computational  effort  required  to  evaluate 
it  is  completely  negligible. 

An  individual  molecule  of  species  i  will  suffer  collisions 
with  molecules  of  species  j  with  a  frequency  given  by 


1 


n 


a .  .c 
ID  r 


(16) 


where  n^^  is  the  number  density  of  species  j  and  is  the 

average  product  of  cross  section  times  relative  velocity  for  the 
two  species,  obtained  by  integrating  over  the  distribution  func¬ 
tion  given  in  Eg.  (1) .  When  this  operation  is  performed,  the 
result  is 
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V? 

1 


(17) 


n. 


r(2  -  (0)  a?.-  , 


ID 


where  T  denotes  the  gamma  function. 

The  total  collision  frequency  for  an  individual  molecule  of 
species  i,  v^,  is  obtained  by  summing  Eq.  (15)  over  all  species, 

i.e. 


1  Z— (  1 

j=i 


(18) 


and  the  mean  free  path,  for  molecules  of  species  i  is  given  by 


/SR.T^ 

A.  =  5T/v.  =  / — - — /v.  ,  (19) 

1  1  Trm^  j  X  ' 

where  cT  is  the  mean  molecular  speed  for  species  i  molecules. 

The  mean  free  path  for  the  gas  mixture,  X^,  is  then  defined  as 
the  number  density  weighted  average  of  the  X^^  via 


n .  X 


Ei<X)  i 
— 

<30 


(20) 


i=l 


where  n  is  the  total  number  density: 


^oo  ”  ^  ,  ^ioo 


i=l 

A  useful  velocity  scale  is  given  by  v  ,  defined  by 

o 


r 


(21) 


where  m  is  the  reference  mean  molecular  weight,  i.e. 


m 


E 

■i 


n , 


m. 
“  1 


n. 


(23) 


Vg  is  the  most  probable  molecular  speed  for  molecules  of  the 
mean  molecular  weight  at  the  reference  temperature. 


2 . 4  Internal  Energy  Model 

The  current  state  of  modeling  for  internal  energy  effects 

in  Monte  Carlo  flow  field  simulations  is  the  phenomenological 

(4) 

model  of  Borgnakke  and  Larsen.  In  this  model,  transfer  of 

energy  between  internal  and  translational  modes  is  allowed, 
but  it  is  necessary  to  assume  that  each  species  has  a  fixed 
number  of  internal  degrees  of  freedom,  This  is  equiva¬ 

lent  to  assuming  a  constant  specific  heat,  C  . ,  for  each  species 
which  can  be  related  to  the  number  of  internal  degrees  of  free¬ 
dom  via 


C  .m. 

C.  =  2  - 

^0 


(24) 


Alternatively,  can  be  related  to  the  ratio  of  specific  heats 
for  species  i,  by  the  relation 


(25) 


The  interchange  of  internal  and  translational  energy  will  be 
discussed  in  Section  4,  and  the  selection  of  initial  conditions 
will  be  considered  in  Section  7. 
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3.  INTERNAL  REPRESENTATION 


3.1  State  Vector 

Each  simulated  molecule  in  the  direct  simulation  Monte  Carlo 
method  is  represented  by  a  state  vector  which  comprises  all  of 
the  information  the  code  has  with  regard  to  that  particular 
molecule.  The  state  vector  has: 

•  Position  element (s)  defining  the  location  of  the 
molecule  in  the  coordinate  system  being  used.  For 
axisymmetric  simulations,  this  is  a  radial  and  an 
axial  element. 

•  Three  velocity  elements.  A  molecular  collision  is 
always  considered  as  a  three  dimensional  event, 
regardless  of  the  overall  dimensionality  of  the 
problem.  For  spatially  one  dimensional  problems 

it  is  possible  to  store  only  two  pieces  of  velocity 
information  and  compute  the  required  three  velocity 
components  as  needed  for  collision  sampling.  This 
is  an  example  of  the  frequent  tradeoff  which  must 
be  made  between  storage  and  computing  requirements. 

•  A  value  for  the  internal  energy  of  the  molecule. 

Note  that  the  basic  model  does  not  discriminate 
between  internal  modes  for  a  particular  species. 

This  can  be  done,  if  desired,  by  introducing 
separate  species  for  the  separate  modes. 

•  An  indicator  determining  the  molecular  species. 

This  indicator  in  turn  implies  all  of  the  proper¬ 
ties  associated  with  that  species  (molecular  weight, 
n;imber  of  internal  degrees  of  freedom,  name,  etc.). 
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•  An  indicator  giving  the  cell  in  which  the  molecule 
currently  resides.  It  is  possible  to  avoid  allo¬ 
cating  this  particular  storage  element,  but  it 
usually  results  in  enough  computational  simplifica¬ 
tion  to  justify  its  use. 


3.2  Reduction  to  a  Reasonable  Number  of  Simulated  Molecules 

It  is  clearly  impossible  to  run  a  computer  simulation  with 
anywhere  near  the  same  number  of  molecules  that  exist  in  the 
actual  flow  problem.  The  adjustment  that  is  made  to  make  the 
simulation  possible  is  to  artificially  increase  the  cross  sec¬ 
tion,  and  decrease  the  number  density,  by  a  large  factor.  It 
is  the  product  of  number  density  and  cross  section  which  deter¬ 
mines  the  collision  frequency  for  a  given  molecule,  and  it  is 
the  collision  frequency  which  must  be  correctly  simulated  if  the 
correspondence  between  the  real  and  simulated  flows  is  to  be 
correct.  This  is  an  essential  feature  of  the  direct  simulation 
method  which  has  not  always  been  adequately  emphasized.  It 
means  that  the  internal  scaling  factors  do  not  proceed  on  a 
strictly  dimensional  basis.  For  example,  the  scaling  factor 
for  cross  sections  is  not  the  square  of  the  scaling  factor  for 
lengths. 


3.3  Internal  Scales 

Many  problems  are  more  reasonably  handled  if  the  internal 
calculations  are  carried  out  with  scaled  or  dimensionless  values. 
This  avoids  possible  problems  such  as  numerical  overflow  which 
can  cause  an  execution  time  error.  Such  errors  can  be  particu¬ 
larly  insidious  and  difficult  to  locate  in  a  code  whose  very 
essence  involves  the  random  combination  of  numbers.  Furthermore, 
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examination  of  scaled  values  makes  the  detection  of  erroneous 
values  easier  while  debugging  codes  since  large  exponents  are 
usually  indicative  of  an  error  when  the  variables  are  internally 
scaled.  At  SSI,  at  least,  the  output  is  produced  in  physically 
meaningful  dimensional  form.  Hence,  the  scaling  that  is  discussed 
here  is  irrelevant  (or  nearly  so)  to  the  interpretation  of  code 
output;  it  is  strictly  a  matter  of  the  internal  representation. 

The  obvious  choices  for  length  and  velocity  scales  are 

and  V  as  defined  in  Section  2,  which  are  used  to  nondimension- 
s 

alize  the  position  and  velocity  elements  of  the  state  vector. 

There  is  no  need  to  provide  further  nondimensionalization  of 

mass  beyond  representing  them  in  atomic  mass  units,  so  none  is 

2 

provided.  Hence,  the  scaling  factor  for  energy  is  just  v  ,  which 

5 

is  used  to  nondimensionalize  the  internal  energy  element  of  the 
state  vector. 

Number  densities  are  scaled  with  respect  to  the  far  field 
reference  number  density,  n^,  which  leaves  only  the  cross  sec¬ 
tion  scaling  factor  to  be  determined.  This  factor  follows  from 
the  condition  of  flow  similarity,  which  requires  that  the  prob¬ 
ability  of  a  molecule  suffering  a  collision  in  traveling  a  given 
path  lem_  _h  be  accurately  simulated.  This  dimensionless  proba¬ 
bility  can  be  expressed  as  the  product  of  a  cross  section  times 
a  nvimber  density  times  a  path  length  (at  least  for  small  enough 
path  lengths) ,  and  it  is  required  that  this  product  be  the  same 
for  dimensional  and  scaled  representations.  This  implies  that 
the  product  of  the  scaling  factors  for  these  three  quantities  be 
unity  and,  therefore,  the  cross  section  scaling  factor  is 
1/ (n  X  ).  The  internal  scaling  factors  used  for  the  SSI  Monte 
Carlo  codes  are  summarized  in  Table  1. 
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TABLE  1  -  SCALING  FACTORS  USED  FOR  THE  INTERNAL 

REPRESENTATION  OF  QUANTITIES  IN  THE  SSI 
DIRECT  SIMULATION  MONTE  CARLO  CODES. 

ALL  VARIABLES  ARE  DEFINED  IN  SECTION  2. 


PROPERTY 

SCALING  FACTOR 

Length 

^00 

Velocity 

^S 

Time 

^co/^S 

Number  Density 

Mass 

a.m.u. 

Energy 

(a.m.u. )Vg 

Cross  Section 

1/ 

3 . 4  Weighting  Factors 

Weighting  factors  are  a  crucial  element  of  a  successful 
Monte  Carlo  simulation,  allowing  trace  species  to  be  described 
with  reasonable  statistics.  The  weighting  factor  is  the  number 
of  "real”  molecules  that  corresponds  to  each  "simulated”  molecule. 
A  "simulated"  molecule  corresponds  to  one  molecule's  worth  of 
storage  (one  state  vector)  allocated  in  the  computer,  and  the 
weighting  factor  is  its  statistical  weight.  So,  for  example, 
the  total  number  density  in  a  cell  might  be  represented 

V 

^cell  “  ~v~ 

i=l 


(26) 
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where  indicates  the  number  of  simulated  molecules  of  species 
i  in  the  cell,  is  the  weighting  factor  for  that  species  in 
that  cell,  V  is  the  cell  volume  and  p  is  the  number  of  species. 

The  product  that  appears  in  Eq.  (26)  is  termed  the  number 

of  "real"  molecules  of  species  i  in  the  cell.  Note  that 
as  calculated  by  Eq.  (26)  is  a  scaled  value;  it  would  have  to 
be  multiplied  by  n^,  as  shown  in  Table  1,  to  become  a  dimensional 
evaluation  of  the  number  density. 

The  weighting  factors  used  in  the  SSI  codes  are  dependent 
on  cell  and  species.  Hence,  flowfields  where  a  given  species 
is  much  more  dominant  in  one  portion  of  the  solution  region  than 
another  can  be  accurately  represented.  It  is  possible,  of  course, 
to  make  weighting  factors  functions  of  other  variables,  such  as 
velocity,  for  specialized  purposes. 

A  critical  error  that  can  occur  in  Monte  Carlo  codes  is  to 
have  the  number  of  simulated  molecules  exceed  the  dimensioned 
limit  of  the  code.  On  the  other  hand,  it  is  generally  desirable 
to  have  as  many  molecules  as  is  feasible  to  obtain  good  statis¬ 
tics.  Resolution  of  these  conflicting  considerations  is  compli¬ 
cated  by  lack  of  a  priori  knowledge  of  what  the  species  number 
densities  will  be  as  a  function  of  space  and  time.  The  way  the 
resolution  is  achieved  in  the  SSI  codes  is  by  a  dynamic  adjust¬ 
ment  of  the  weighting  factors,  as  required.  This  keeps  the 
number  of  simulated  molecules  more  or  less  constant  while 
allowing  the  number  of  real  molecules  to  adjust  as  the  solution 
evolves.  The  introduction  of  weighting  factors,  with  the  abil¬ 
ity  to  adjust  them  as  the  solution  demands,  is  an  important 
feature  of  a  successful  Monte  Carlo  description. 
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4.  COLLISION  MECHANICS 


4 . 1  Relations  for  Elastic  Collisions 

The  purpose  of  this  section  is  to  present  relations 
appropriate  to  the  simulation  of  a  collision  in  a  Monte  Carlo 
flowfield  code.  (The  question  of  how  molecules  are  selected  for 
collisions,  which  is  crucial  to  the  proper  simulation  of  colli¬ 
sion  frequency,  will  be  taken  up  in  Section  6.)  Conservation  of 
momentum  implies  that  the  center-of-mass  velocity  of  the  colli¬ 
sion  pair  is  unchanged  by  the  collision;  and  conservation  of 
energy  then  implies  that  the  magnitude  of  the  relative  velocity 

between  the  collision  partners  is  also  unchanged  by  the  colli- 
15) 

Sion.  '  Since  the  collision  is  treated  as  a  statistical  event, 
all  that  remains  is  to  select  the  direction  of  the  post-collision 
relative  velocity  vector  from  the  correct  distribution.  As 
mentioned  in  Section  2.2,  collisions  ii  the  VHS  model  are 
treated  as  hard  sphere  collisions  when  they  occur  (though  they 
do  not  occur  with  the  same  velocity  dependence  as  do  hard  sphere 
collisions) .  Hence,  as  far  as  the  collision  mechanics  is  con¬ 
cerned,  the  model  is  a  hard  sphere  model.  For  hard  sphere 
molecules,  all  directions  for  the  post-collision  relative 
velocity  vector  are  equally  likely.  This  is  the  chief  computa¬ 
tional  simplicity  of  the  VHS  model. 

Let  the  two  molecules  be  identified  by  subscripts  1  and  2, 
with  m  and  v  denoting  their  masses  and  velocities.  If  i  and  f 
indicate  initial  and  final  states,  then  the  relations  for  the 
collision  can  be  summarized  via: 
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(27) 


^cm 

“l^li  +  ">2^21 

i"l  ">2 

(27) 

^r 

l^li  -  ^2i|  ' 

(28) 

cos  (0) 

=  1  -  2B 

(29) 

sin  (9 ) 

=  (9)  , 

(30) 

(|)  = 

27rB 

(31) 

II 

M-l 

v^^cos(9),  sin  (9)cos  ((j>)  ,  sin  (0)  sin  ((j))J 

f 

(32) 

^cm  +  m2  ^rf  ' 

(33) 

^2f  = 

m  ^ 

V  —  -  V  ^  . 

cm  mj^  +  m2  rf 

(34) 

In  these  relations,  and  throughout  this  report,  3  indicates 
a  random  variable  which  is  evenly  distributed  on  the  interval 
zero  to  one.  Each  time  that  B  appears  a  different  evaluation  of 
the  random  variable  is  implied.  Note  that  the  expression  for  the 
post-collision  relative  velocity  vector  (Eq.  (32))  is  not  coordi¬ 
nate  system  specific.  The  indicated  vector  components  can  apply 
to  any  locally  orthogonal  coordinate  system,  since  the  direction 
implied  is  random.  The  convenient  coordinate  system  to  use,  of 
course,  is  the  coordinate  system  used  to  define  the  velocity 
elements  of  the  state  vector.  For  axisymmetric  simulations  this 
will  be  radial,  azimuthal  and  axial  velocity  components. 
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4.2  Simulation  of  Inelastic  Collisions 


(4 ) 

The  SSI  codes  use  the  Borgnakke  and  Larsen  phenomenolog¬ 
ical  model  for  transfer  of  energy  between  internal  and  transla¬ 
tional  modes.  In  this  model,  a  collision  is  assumed  to  be  either 
perfectly  elastic  or  perfectly  inelastic,  via  a  user  specified 
probability.  A  perfectly  inelastic  collision  is  achieved  by 
summing  the  total  pre-collision  energy  (internal  energy  of  both 
molecules  plus  the  translational  energy  of  their  relative  motion, 
Eq.  (4))  and  then  assigning  post-collision  values  from  the 
equilibrium  distribution  for  collisions  with  that  total  amount 
of  energy,  taking  into  account  the  number  of  internal  degrees  of 
freedom  in  the  two  molecules.  Note  that  this  model  has  the 
ability  to  relax  from  a  nonequilibrium  to  an  equilibrium  state 
via  an  effective  collision  number.  The  ability  to  exchange 
internal  energy  in  such  a  manner  comprises  a  significant  increase 
in  capability  for  Monte  Carlo  codes  beyond  the  previous  models 
where  molecules  had  no  internal  energy.  It  is  this  capability 
which  enables  the  codes  to  realistically  predict  the  macroscopic 
effects  of  polyatomic  gas  flow. 

Let  and  ^2  the  number  of  internal  degrees  of  freedom 

of  the  two  molecules  in  an  inelastic  collision,  and  E  be  the 

s 

total  collision  energy  defined  by 


=  E  .  + 
Cl 


Eli  +  E2i 


(35) 


where  E^^  is  the  initial  translational  collision  energy  defined 
by  Eq.  (4),  and  E^^  and  pre-collision  internal  ener¬ 

gies  of  the  two  molecules.  If  C  is  defined  by 

5  =  ,  (36) 


where  E^^  is  the  post-collision  translational  energy,  then  C  is 
selected  according  to  the  distribution 


f  (C) 


(37) 


where 


and 


A  = 

r<c>  - 

i-“r<c>  - 

[l  -  (U  J 

L<^>  -  ij 

b  = 

<?>  - 1 

<?> 

II 

+ 

^2)/2 

The 

sampling  of 

Eq.  (37)  is 

<C>-1 


(38) 


(39) 


(40) 


used  frequently  in  Monte  Carlo  flowfield  codes,  the  acceptence- 
rejection  method.  Equation  (37)  has  been  normalized  so  that  the 
maximiim  value  of  the  function  is  unity,  and  the  parameter  of  the 
distribution,  5,  varies  from  zero  to  one.  The  sampling  is  done  as 
follows: 


Choose  a  random  value  of  I.e.,  set  C  equal  to 
6,  a  random  variable. 

Evaluate  the  distribution  function  for  this  value 
of  5/  and  call  it 

Get  a  second  random  variable,  and  check  to  see  if 
it  is  greater  or  less  than  f^gg^*  If  if  is  greater 
than  f^gg^/  9°  back  to  the  first  step  and  repeat 
the  process.  If  it  is  less  than  f^gg^/  then  keep 
the  value  of  ?  obtained  in  the  first  step. 


Note  that  the  probability  that  any  original  value  of  C  will  be 
kept  is  proportional  to  f(C).  This  is  an  extremely  general  tech¬ 
nique  and,  as  such,  is  very  powerful.  It  is  not  always  efficient, 
however,  and  direct  sampling  of  distributions  is  usually  to  be 
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preferred  if  it  can  be  accomplished.  In  this  case,  if  one  mole¬ 
cule  is  monatomic  and  the  other  is  diatomic,  i.e.  if  <c;>  =  1, 
then  a  direct  sampling  of  the  above  distribution  is  achieved  via 

Once  C/  and  therefore  the  post-collision  translation  energy, 
is  determined,  the  magnitude  of  the  post-collision  relative 
velocity  is  defined  via 

■  ''r  =  •  '‘'2> 

For  inelastic  collisions,  this  relation  takes  the  place  of  Eq. 

(28)  in  the  determination  of  the  post-collision  velocity  elements 
of  the  state  vectors. 

The  remainder  of  the  collision  energy  must  be  divided  up 
between  the  internal  modes  of  the  two  molecules.  If  one  of  the 
molecules  is  monatomic  (i.e.,  has  zero  internal  degrees  of  free¬ 
dom),  then  all  of  the  internal  energy  goes  to  the  other  by  default. 
Otherwise,  if  x  is  defined  by 

where  is  the  post-collision  internal  energy  of  the  first 

molecule,  then  x  is  just  the  fraction  of  the  total  post-collision 
internal  energy  that  ends  up  in  the  first  molecule,  x  is  distri¬ 
buted  with  a  probability  proportional  to  f(x),  given  by 

f(x)  =  B  x'^d-x)^  ,  (44) 

where 

B  =  [(<?>-  2)  ,  (45) 

c  =  Qj^/2  -  1  ,  (46) 
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and 


(47) 


d  =  C2/2  -  1 

Equation  (44)  can  be  sampled  via  the  acceptance-rejection  method 
to  determine  the  allocation  of  the  internal  energy  for  the  general 
case.  For  the  special  case  of  both  molecules  being  diatomic,  Eq. 
(44)  becomes  singular,  with  the  limit  being  the  trivial  case  that 
all  distributions  of  internal  energy  are  equally  likely,  i.e. 

X  =  6  ,  (;^  =  =  2)  .  (48) 

For  the  special  case  of  just  one  of  the  molecules  being  diatomic 
(molecule  2,  for  instance),  then  the  distribution  for  x  can  be 
sampled  directly  via 
(2/C,) 

X  =  6  ,  (C2  =  2)  ,  (49) 

with  the  obvious  reciprocal  relation  applying  for  =  2.  The 
SSI  codes  recognize  these  special  cases  so  that  the  sampling  can 
be  expedited  when  possible,  while  retaining  the  full  generality 
of  Eqs.  (37)  and  (44)  when  required. 

4 . 3  Collisions  for  Molecules  with  Distinct  Weighting  Factors 

There  is  an  obvious  problem  when  considering  a  collision 
between  two  simulated  molecules  with  distinct  weighting  factors, 
since  they  represent  a  different  number  of  real  molecules.  If 
Wy  and  Wy  represent  the  weighting  factors  for  the  two  molecules, 
with  being  greater  than  ,  then  the  collision  is  always 
counted  as  "events”.  This  ii  accomplished  by  always  assign- 
ing  post-collision  velocity  and  energy  components  to  the  state 
vector  of  the  molecule  with  the  smaller  weighting  factor,  but 
only  changing  the  components  of  the  molecule  with  the  greater 
weighting  factor  some  of  the  time.  The  probability  that  the 
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molecule  with  the  greater  weighting  factor  will  have  its  compon¬ 
ents  changed  is  simply  W  /W  .  Statistically,  this  means  that 
for  a  large  niamber  of  simulated  collisions,  each  such  simulated 
collisioii  will  average  out  to  W  real  collisions  for  each  species, 
even  though  their  weighting  factors  differ.  It  should  be  noted 
that  this  does  violate  conservation  of  momentum  and  energy  on  an 
individual  collision  basis,  but  these  quantities  are  conserved  in 
the  aggregate  over  a  large  number  of  collisions. 


4 . 4  Reactive  Collisions 

A  realistic  simulation  of  chemical  reactions  is  a  crucial 
element  of  many  problems.  If  a  bimolecular  reaction  of  the  form 

A  +  B  ->  C  +  D  (50) 

has  an  Arrhenius  rate  constant  of  the  form 


k  =  A^T^  exp(-E^/RQT) 


(51) 


then  it  is  possible  to  define  a  reactive  cross  section  as  a  func¬ 
tion  of  translational  collision  energy  such  that  the  above  rate 
constant  is  implied  for  a  gas  in  translational  equilibrium.  (In 
Eq.  (51),  T  is  temperature,  E,  is  the  activation  energy,  n  is  a 
dimensionless  exponent  and  A^  is  a  prefactor.)  The  product  of 
reactive  cross  section,  a*,  times  relative  velocity  can  be 
expressed  ^ 


V  a* 
r 


(1  +  6^j)  VT"  A^ 
2Rq  r(n  +  3/2) 


(52) 


where  is  unity  for  like  reactants  and  zero  for  unlike  reac¬ 

tants.  The  probability  of  a  collision  resulting  in  a  reaction  is 
given  by  the  ratio  of  the  reactive  to  the  gas  kinetic  cross  sec¬ 
tion  at  the  existing  relative  velocity  between  the  molecules. 
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In  the  SSI  Monte  Carlo  codes,  the  following  procedure  is 
used  to  determine  which,  if  any,  reaction  will  occur: 

•  All  possible  reactions  are  identified  for  the  pair 
of  molecules  which  are  selected  to  experience  a 
collision.  If  there  are  no  such  reactions, 

the  procedure  is  bypassed. 

•  The  probability  of  each  allowed  reaction  is  calcu¬ 
lated  by  computing  the  reactive  to  gas  kinetic 
cross  section  ratio. 

•  If  somehow  the  reactive  cross  section  (s)  total 
to  a  greater  value  than  the  gas  kinetic  cross 
section  (which  will  rarely  be  the  case)  then 
the  probabilities  are  normalized  by  their  sum. 

•  The  collision  is  taken  to  be  reactive  with  a  prob¬ 
ability  equal  to  the  sum  of  the  individual  reaction 
probabilities . 

•  If  the  collision  is  reactive,  the  reaction  that 
occurs  is  selected  in  accordance  with  its 
probability. 

•  If  no  reaction  is  decided  upon,  then  the  collision 
is  either  elastic  or  inelastic  according  to  the 
user  specified  probability. 

Equation  (52)  becomes  singular  for  n  <  -3/2,  and  this  pro¬ 
cedure  cannot  handle  that  case.  The  major  limitation  of  this 
procedure  is  that  it  ignores  the  effect  of  internal  energy  in 
determining  whether  or  not  a  reaction  can  occur,  assuming  that 
the  entire  activation  energy  barrier  must  be  overcome  through 
translational  energy.  If  a  specified  fraction  of  the  molecular 
internal  energy  is  allowed  to  contribute  to  overcoming  this 
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barrier,  then  the  restriction  on  n  is  relaxed  somewhat  (see 
Ref.  3),  but  this  feature  is  not  presently  implemented  in  the 
SSI  codes. 

4 . 5  Post-Reaction  State  Definition 

The  post-reaction  state  is  simply  determined  by  adding  the 
heat  of  reaction  to  the  total  collision  energy  (translational 
plus  internal)  and  then  computing  reactant  state  vectors  for  the 
products  based  on  this  total  energy  and  the  reactants  center-of- 
mass  velocity,  as  for  inelastic  collisions.  Although  this  is 
clearly  an  oversimplification,  it  is  quite  consistent  with  the 
phenomenological  nature  of  the  model. 

The  position  state  vector  components  for  the  products  are 
randomly  selected  from  the  position  state  vector  components  of 
the  reactants,  which  are  not  the  same  (see  Section  6) .  The 
major  additional  complication  of  chemical  reactions  is  that  of 
distinct  weighting  factors.  Since  the  reaction  is  "events" 
(see  Section  4.3),  it  only  destroys  the  reactant  with  the  greater 
weighting  factor  with  a  probability  of  If  Wp  is  a  product 

weighting  factor,  it  is  necessary  to  produce  W^/W  simulated 

Jj  p 

molecules  of  that  product.  In  general  this  is  not  an  integer 
quantity,  so  it  is  necessary  to  interpret  the  ratio  statistic¬ 
ally  so  that  the  expectation  value  is  proper.  That  is,  sometimes 
the  next  lower  integer  is  selected  and  sometimes  the  next  higher 
one,  with  a  probability  that  reflects  how  close  each  integer 
value  is  to  the  desired  fractional  quantity.  (See  the  discussion 
of  molecular  cloning  in  Section  5.2.)  Of  course,  the  weighting 
factors  for  the  two  products  will  in  general  be  different, 
resulting  in  a  different  number  of  simulated  molecules  being  pro¬ 
duced  for  the  two  products.  Sometimes  this  process  could  result 


in  a  very  large  production  of  simulated  molecules  for  a  product 
with  a  small  weighting  factor.  In  order  to  prevent  overflow  of 
code  dimensions  it  is  necessary  for  the  code  to  sense  when  this 
is  happening  and  automatically  increase  the  product  weighting 
factor  to  prevent  it. 


4.6  Dissociative  Reactions 

The  situation  for  dissociative  reactions  is  somewhat  compli¬ 
cated  by  the  presence  of  three  rather  than  two  products.  With  a 
little  manipulation,  however,  it  is  possible  to  use  the  previous 
relations  for  this  case  as  well.  Let  M^,  M2  and  represent  the 
three  product  molecules  from  a  dissociative  reaction,  with  a 
known  center-of-mass  velocity  and  total  energy.  The  procedure 
for  defining  the  post-reaction  state  is  to  first  define  an  arti¬ 
ficial  complex  comprised  of  the  (M2»M2)  pair.  (There  is  no 
implication  that  there  actually  is  any  such  collision  complex.) 
The  complex  is  assigned  a  number  of  internal  degrees  of  freedom 
equal  to  5^,  given  by 

?c  “  ?2  ^3  2  (2  -  to)  ,  (53) 

where  the  last  term  represents  the  contribution  of  the  relative 
translation  between  M2  and  M^. 

The  fact  that  this  term  is  not  simply  three,  as  it  would 
seem  it  should  be,  merits  some  explanation.  It  is  due  to  the 
fact  that  these  molecules  are  not  random  samples  from  the  gas  but 
rather  special  molecules  owing  to  their  being  created  in  a  reac¬ 
tion.  This  point  can  perhaps  best  be  seen  by  considering  micro¬ 
scopic  reversibility,  where  the  inverse  reaction  is  a  three  body 
recombination.  For  this  reverse  process,  molecules  participating 
in  it  are  not  all  equally  probable,  since  those  with  greater 
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relative  velocities  are  more  likely  to  collide.  Hence,  the  term 
does  take  on  the  value  three  for  the  special  case  of  oo  equal  to 
1/2,  which  is  precisely  the  case  of  collision  frequency  being 
independent  of  relative  velocity.  Examination  of  Eqs.  (37)  and 
(44)  demonstrates  that  translational  energy  in  collisions  behaves 
like  another  source  of  internal  energy  with  2(2  -  u)  degrees  of 
freedom. 

With  the  number  of  degrees  of  freedom  defined,  the  separa¬ 
tion  of  from  the  (M2,M2)  complex  is  treated  as  an  inelastic 
collision.  The  resulting  velocity  for  the  complex  is  then 
treated  as  the  center-of-mass  velocity  for  M2  and  M^/  and  the 
internal  energy  assigned  to  the  complex  becomes  the  total  energy 
for  the  pair.  Using  these  values,  M2  and  M^  are  then  separated 
in  another  application  of  the  rules  for  inelastic  collisions. 


5.  MOLECULAR  TRANSLATIONS 


5.1  Translation  in  an  Axi symmetric  Coordinate  System 

As  discussed  in  Section  1,  an  essential  element  of  the  direct 
simulation  Monte  Carlo  method  is  the  periodic  advancement  of  sim¬ 
ulated  molecules  along  their  trajectories.  Formally,  this  is 
accomplished  by  updating  the  position  and  velocity  elements  of 
the  state  vector.  For  a  cartesian  coordinate  system  this  is  a 
trivial  process,  but  the  relations  are  slightly  more  complicated 
for  the  often  used  axisymmetric  coordinate  system.  Let  v^q, 
v^Q,  and  v^Q  represent  the  initial  radial,  azimuthal  and  axial 
velocity  components  of  a  molecule,  with  r^  and  z^  representing 
its  initial  radial  and  axial  position.  Additionally,  let  <J)q 
represent  the  initial  azimuthal  angle  for  the  molecule.  This  is 
included  here  just  for  demonstration  purposes;  it  will  not  gen¬ 
erally  be  known,  nor,  as  will  be  shown,  will  it  be  needed. 

The  initial  position  and  velocity  of  the  molecule  can  then 
be  referenced  to  a  standard  cartesian  coordinate  system,  yielding 

^0  =  [^rO  cos((Dq)  -  v^Q  sin(<})Q)ji 

+  [vj.0  sin(<p^)  v^Q  cos(4>o)]j 


=  V 


xO 


1  +  V 


yO 


j  +  V  -  k 
zO 


(54) 
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and 


'O  =  [""o  cos{(j)Q)ji  +  jr^  sin(<i)Q)Jj  +  Zq  k 

^ 

=  Xq  i  +  Yq  j  +  Zq  k 


(55) 


After  a  time  increment  Dt,  the  position  vector  of  the  molecule 
will  be 

=  Dt)  i  +  (y^  +  Dt)  j  +  (z^  +  v^q  Dt)  k 


1  '^0  '  "xO  ^  ■  '^0  "yO 


=  X,  i  +  y,  j  +  z,  k 


(56) 


and,  in  the  absence  of  perturbing  forces,  the  velocity  vector  will 
remain  unchanged  in  the  cartesian  coordinate  system.  It  will 
change  in  the  axisymmetric  coordinate  system  since  the  basis 
vectors  are  a  function  of  position  and  the  molecule  has  moved. 

The  new  radial  position  is 


■1  = 

-  V 


<■^0  '^rO  +  '''♦O 


and  the  new  axial  position  is 


^1  "  ^0  '^zO 


(57) 


(58) 


The  radial  velocity  component  in  the  coordinate  system  appropriate 
to  the  new  molecular  position  can  be  determined  by  noting 


''rl  =  <''0>  • 


-  ['^ 


ro'^^O 


(59) 
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and  it  can  similarly  be  shown  that 


and 


zl 


=  V 


zO 


(60) 


(61) 


Equations  (57)  -  (61)  give  the  updated  position  and  velocity 
elements  of  the  state  vector  for  a  translation  in  an  axi symmetric 
coordinate  system.  Note  that  these  relations  are  indeed  independ¬ 
ent  of  (J)q.  a  similar  procedure  applies  for  molecular  transla¬ 
tion  in  other  coordinate  systems. 

5.2  Molecular  Cloning 

When  a  simulated  molecule  is  translated  from  one  cell  to 
another,  the  weighting  factor  for  that  species  will  generally 
be  different  in  the  new  cell.  Since  it  is  the  number  of  real 
molecules  rather  than  the  number  of  simulated  molecules  which 
must  be  preserved  when  crossing  cell  boundaries  (statistically, 
at  least) ,  it  is  necessary  to  correct  for  the  distinct  weighting 
factors  (see  Section  3.4). 

If  the  weighting  factor  before  translation  is  Wq,  then  the 
simulated  molecule  represents  that  many  real  molecules.  If  the 
weighting  factor  in  the  new  cell  is  W^,  then  Wq/W^^  simulated 
molecules  would  be  required  to  represent  the  same  number  of  real 
molecules  in  the  new  cell.  If  this  ratio  were  a  whole  integer, 
then  this  could  be  accomplished  by  introducing  that  many  "clones" 
of  the  molecule  in  the  new  cell.  That  is,  that  many  simulated 
molecules  would  be  placed  in  the  new  cell,  all  with  the  same 
state  vector. 


When  the  niunber  Wq/W^^  is  not  an  integer  (the  usual  case,  of 
course),  then  the  cloning  must  be  done  on  a  statistical  basis. 

So,  for  instance,  if  Wq/W^  were  equal  to  2.7,  then  30%  of  the 
time  two  clones  would  be  produced  and  70%  of  the  time  three 
would  be  produced.  Note  that  the  ratio  may  be  less  than  unity, 
and  the  molecule  may  not  be  introduced  into  the  new  cell  at  all. 
(In  which  case  the  molecule  is  removed  from  the  simulation.) 

Cloning  is  a  necessary  evil  inherent  in  a  system  with  spa¬ 
tially  dependent  weighting  factors.  It  enables  such  a  system  to 
maintain  the  statistically  correct  flux  of  mass  and  momentum 
across  cell  boundaries,  but  it  misrepresents  the  flux  of  random¬ 
ized  or  thermal  energy.  This  can  be  seen  by  an  extreme  case 
where  a  very  large  number  of  clones  is  produced  when  a  simulated 
molecule  crosses  a  cell  boundary.  The  resulting  molecules  in 
the  new  cell  have  the  correct  mass  and  momentum  flux,  but  since 
they  all  have  precisely  the  same  velocity  they  have  a  null  rela¬ 
tive  velocity,  and  therefore  a  zero  temperature.  If  the  weighting 
factors  are  not  too  different  between  adjaceiit  cells,  then 
the  errors  introduced  by  this  process  are  acceptably  small. 
However,  it  does  mean  that  one  cannot  arbitrarily  improve  statis¬ 
tics  in  one  portion  of  the  solution  region  by  selectively  reduc¬ 
ing  the  weighting  factors  there.  This  was  a  difficulty  which  was 
encountered  in  the  early  stages  of  the  direct  simulation  Monte 
Carlo  method  while  trying  to  improve  statistics  along  the  axis  of 
axisymmetric  simulations,  since  the  cell  volumes  (and  therefore 
the  sample  sizes)  tend  to  be  smallest  on  the  axis. 

As  was  the  case  for  simulated  molecules  produced  via  chemical 
reactions,  it  is  possible  for  the  weighting  factors  between  suc¬ 
cessive  cells  to  be  so  different  that  a  prohibitively  large  number 
of  simulated  molecules  would  be  required  to  produce  the  same 
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number  of  real  molecules.  This  is  most  frequently  the  case  when 
a  new  species  is  being  introduced,  since  before  the  species  gets 
there  the  weighting  factor  is  initialized  to  a  very  small  number. 
The  codes  sense  when  a  disproportionate  number  of  simulated  mole¬ 
cules  are  being  produced  for  a  given  species  and  cell,  and  adjust 
the  weighting  factor  automatically.  As  the  weighting  factor  is 
increased,  a  proportionate  fraction  of  molecules  of  that  species 
and  cell  are  removed  from  the  simulation  in  order  to  keep  the 
number  of  real  molecules  properly  represented.  This  process 
enables  the  weighting  factors  to  seek  their  own  proper  level 
without  a  priori  knowledge  of  the  solution.  (Periodically, 
the  cells  are  examined  to  determine  if  a  certain  species  has 
been  underrepresented  in  terms  of  number  of  simulated  molecules. 
If  it  is  found  to  be  the  case,  then  the  weighting  factor  is 
decreased,  allowing  weighting  factors  to  float  downwards  as  well 
as  upwards.  It  is  the  danger  of  weighting  factors  being  too 
small,  causing  an  overflow  of  code  dimensions,  which  is  most 
critical,  however.) 
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6.  COLLISION  SAMPLING  IN  A  MULTI -COMPONENT  VHS  GAS 

6.1  General  Considerations  and  Approach 

The  two  general  considerations  in  the  sampling  of  collisions 
are,  as  usual,  accuracy  and  efficiency  of  the  simulation.  As 
far  as  accuracy  is  concerned,  it  is  crucial  that  the  method  in 
which  molecules  are  selected  for  collisions  be  proper.  It  is 
imperative  that  the  correct  collision  frequencies  be  simulated 
between  various  species,  and,  in  fact,  between  the  different 
portions  of  the  velocity  phase  space  for  the  various  species. 
Furthermore,  this  frequency  of  simulated  collisions  must  remain 
correct  without  any  requirements  put  on  the  velocity  distribution 
function;  it  certainly  must  not  be  assumed  that  there  is  a 
Maxwellian  velocity  distribution. 

As  far  as  efficiency  is  concerned,  it  is  highly  desirable 
to  use  a  method  of  collision  sampling  involving  a  computational 
effort  which  is  proportional  to  the  number  of  simulated  molecules, 
N,  in  a  cell.  Methods  which  are  proportional  to  a  power  of  N 
greater  than  unity  can  become  prohibitively  time  consuming  as 
the  number  of  molecules  is  increased  —  a  limit  which  should  be 
made  as  accessible  as  possible  for  obvious  physical  reasons. 


6 . 2  Collision  Sampling  for  a  Single  Component  Gas 

The  simplified  situation  of  a  simulation  involving  only  one 
species  is  considered  here.  This  problem  is  significant  in  part 
due  to  all  the  attention  it  has  received  and,  as  will  be  seen, 
it  serves  as  an  important  reference  case.  When  there  is  just 
one  species,  then  there  is  just  one  gas  kinetic  cross  section 
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(though  it  is  still,  of  course,  a  function  of  collision  energy), 
just  one  molecular  weight  and  just  one  weighting  factor  for  each 
cell.  In  short,  just  one  of  everything  that  has  a  molecular  sub¬ 
script.  Hence, in  this  subsection  all  such  quantities  will  be 
presented  without  subscripts.  The  most  important  simplification 
of  having  a  single  species  is  that  there  is  just  one  collision 
class,  i.e.,  only  self-collisions  of  the  given  species  with  it¬ 
self  are  possible. 


6.2.1  Collision  Pair  Selection 

As  discussed  in  Section  1,  in  the  direct  simulation  Monte 
Carlo  method  collisions  are  sampled  on  a  cell  by  cell  basis 
until  the  number  of  collisions  simulated  is  appropriate  to  the 
overall  solution  time  step,  At^.  The  only  spatial  requirement 
placed  on  potential  collision  partners  is  that  they  be  within 
the  same  cell.  In  particular,  it  is  not  required  that  they  be 
within  a  molecular  diameter  of  each  other.  (Note  that  if  all 
pairs  of  molecules  were  inspected  to  find  those  that  were  suffi¬ 
ciently  close  to  each  other,  this  would  involve  a  computational 
effort  in  proportion  to  the  square  of  the  number  of  molecules  in 
the  cell.)  The  rationale  for  this  is  that  the  cells  should  be 
taken  small  enough  such  that  macroscopic  properties  can  be 
assumed  constant  across  the  cell.  When  this  is  the  case,  then 
a  molecule  within  the  cell  can  be  considered  typical  of  a  mole¬ 
cule  which  might  exist  anywhere  within  the  cell,  and  molecular 
location  can  be  ignored  when  selecting  potential  collision  pairs. 

Spatial  consideration  aside,  the  probability  of  any  two 
molecules  experiencing  a  collision  is  proportional  to  ac^,,  the 
product  of  their  mutual  cross  section  times  their  relative  velo¬ 
city.  This  probability  is  correctly  simulated  via  an  application 
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of  the  acceptance-rejection  method  if  pairs  of  molecules  are 
selected  at  random,  and  then  kept  or  rejected  as  collision 
partners  with  a  probability  proportional  ac^.  This  is  accomp¬ 
lished  by  keeping  a  maxim\im  value  for  ac^  for  each  cell  (which 
is  updated  if  a  larger  value  is  encountered)  and  computing  the 
ratio  r  defined  by 


r 


°°r 

'°=r'max 


(62) 


A  random  variable,  6,  is  then  generated,  and  the  pair  of  molecules 
is  accepted  as  collision  partners  if  r  is  greater  than  0.  This 
produces  the  proper  relative  collision  probability  without  regard 
to  the  existing  velocity  distribution  function. 


6.2.2  Collision  Time  Counter  for  a  Single  Component  Gas 

The  volumetric  collision  frequency  for  a  single  component 
gas,  V  (collisions  per  unit  volume  per  unit  time),  is  given  by 


where,  as  in  Section  2,  n  represents  the  number  density  of  the 

species  and  ac^  is  the  average  product  of  collision  cross  section 

and  relative  velocity.  At  first  inspection,  it  would  seem  from 

Eq.  (63)  that  a  correct  simulation  of  collision  frequency  would 

require  evaluation  of  ac^,  which  would  mean  that  all  pairs  of 

molecules  in  a  cell  would  have  to  be  considered.  Such  a  proced- 

2 

ure  involves  a  computational  effort  proportional  to  N  and  is  to 
be  avoided,  if  possible,  in  preference  to  a  method  which  is 
simply  proportional  to  N. 
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The  alternative  approach,  introduced  by  Bird,  is  the  time 
counter  approach.  For  each  collision  a  time  counter,  t^,  is 
incremented  by  an  amount  that  depends  on  the  relative  velocity 
of  the  collision.  Collision  sampling  continues  in  a  cell  until 
its  time  counter  has  been  advanced  beyond  the  overall  flow  simu¬ 
lation  time,  at  which  time  the  code  proceeds  to  the  next  cell 
(which  has  its  own  time  counter).  The  time  counter  increment, 
At^  is  given  by 


At. 


Vn^ac, 


(64) 


where  V  is  the  cell  volume  and  n  is  the  species  number  density 
given  by 

n  =  NW/V  ,  (65) 


with  W  being  the  weighting  factor  for  the  species.  (Equation 
(65)  is  just  a  special  case  of  Eq.  (26).)  It  should  be  stressed 
that  Eq.  (64)  applies  for  each  real  collision.  As  is  discussed 
in  Sections  3.4  and  4.3,  each  simulated  collision  corresponds  to 
W  real  collisions,  so  when  a  simulated  collision  occurs  the 
actual  applied  increment  to  t^  is  W  times  the  value  given  by 
Eq.  (64 ) . 

It  is  not  obvious  that  Eq.  (64)  will  lead  to  a  proper  simula¬ 
tion  of  the  overall  collision  frequency,  so  a  demonstration  will 
be  presented.  Let  f^  normalized  distribution  function 

for  relative  velocity  appropriate  to  a  given  cell  at  a  given  time. 
(Note  that  most  problems  solved  by  a  Monte  Carlo  technique  involve 
repeated  runs,  where  the  total  number  of  collision  pairs  can  be 
made  arbitrarily  large.  The  introduction  of  a  distribution  func¬ 
tion,  and  therefore  the  demonstration  of  the  correctness  of  Eq. 

(64)  is  strictly  valid  only  in  the  limit  of  infinitely  many  runs. 
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since  the  essence  of  the  method  is  a  correct  statistical  repre¬ 
sentation,  it  could  not  be  otherwise.)  By  definition. 


/■ 


f  {Cj.)dCj 


=  1 


(66) 


irrespective  of  the  particular  form  for  f^.  Let  f2(c^)  be  the 
normalized  distribution  function  of  relative  velocity  occurring 
in  collisions.  Since  collisions  occur  with  a  probability  pro¬ 
portional  to  ac^,  f2  can  be  constructed  from  f^  via 


ac^ 

f2(Or)  =  — 

cc_ 


(67) 


where  ac  can  now  be  formally  defined  via 

■  f\ 


(Cj.)oc^dc^ 


(68) 


(Note  that  in  Eq.  (68),  as  in  the  rest  of  this  demonstration,  the 
functional  form  of  the  dependence  of  cross  section  on  relative 
velocity  need  never  be  specified.  The  time  counter  represented 
by  Eq.  (64)  is  therefore  not  restricted  to  any  particular  model 
for  the  cross  section.) 

The  average  increment  of  the  time  counter  that  is  applied 
over  many  simulated  collisions  is  therefore  given  by 


■  r 


f2(Cr)Atc(Cr)dCr 


(69) 


If  Eqs.  (64)  and  (67)  are  substituted  into  Eq.  (69),  the  result  is 

2 


At, 


n  Vac 


(70) 
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where  the  normalization  condition  (Eq.  (66))  has  been  utilized. 
The  implication  of  Eq.  (70)  is  that,  on  average,  the  frequency 
of  collisions  in  the  cell  is  Vn  ac^/2.  Since  this  is  simply 
the  product  of  cell  volume  times  the  proper  volumetric  collision 
frequency  (Eq.  (63)),  the  validity  of  the  time  counter  given  in 
Eq.  (64)  has  been  demonstrated. 


6.3  Collision  Class  Sampling  in  Gas  Mixtures 


The  above  procedure  for  a  single  species  gas  can  be  extended 
to  a  multi-component  mixture  via  consideration  of  distinct  colli¬ 
sion  classes.  In  this  approach,  collision  classes  are  defined 
by  the  colliding  pair  identities.  Hence,  if  there  are  p  species 
in  the  simulation  then  there  are  p(p+l)/2  collision  classes, 

which  can  be  identified  by  the  subscripts  of  the  corresponding 

2 

molecular  pair.  (The  number  of  classes  is  not  p  since  the 
order  of  molecule  specification  is  not  taken  to  matter  in 
determining  a  collision  class.  Hence,  the  class  identified  by 
the  subscripts  i,j  is  not  distinct  from  the  class  identified 
by  the  subscripts,  j,i.) 

In  collision  class  sampling,  which  is  the  method  used  by 
Bird, each  collision  class  is  sampled  separately,  and  the 
collision  sampling  in  a  cell  is  not  complete  until  all  classes 
have  been  considered.  Each  collision  class  has  its  own  stored 
value  of  (j. -c^)  ^  ,  and  its  own  separate  time  counter,  t  • 

By  a  comparable  analysis  to  that  presented  above,  it  can  be 
shown  that  the  appropriate  time  counter  increment  in  this  case  is 


At 


Cl] 


(71) 


where,  as  in  Section  4,  6^^  is  the  Kronecker  delta  which  is 
unity  for  i  =  j  and  zero  otherwise.  As  in  the  previous  section. 
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the  above  increment  applies  for  each  real  collision.  A  simulated 

collision  corresponds  to  W-  real  collisions,  where  is  the  les- 

ser  of  VI ^  and  (see  Section  4.3),  so  when  a  simulated  collision 

occurs,  the  applied  increment  to  t  . .  is  W_  times  the  result  of 
^  C13  L 

Eq.  (71). 


6.4  Global  Collision  Sampling  in  a  Gas  Mixture 

Although  the  procedure  described  above  is  quite  reasonable 
for,  say,  a  two  component  mixture,  it  becomes  quite  complicated 
as  the  number  of  species  increases.  For  10  species,  for  instance, 
the  program  must  loop  over  55  distinct  collision  classes  for  each 
cell,  and  storage  must  be  allocated  for  110  quantities  in  each 
cell.  As  the  number  of  species  increases,  the  storage  require¬ 
ment  for  the  collision  sampling  constants  quickly  becomes  greater 
than  the  storage  required  for  the  molecular  state  vectors!  The 
obvious  simplification  is  to  search  for  a  technique  where  colli¬ 
sions  are  simulated  simultaneously  for  all  collision  classes, 
with  each  class  having  its  proper  relative  probability  of  being 
selected.  The  overall  collision  sampling  then  continues  until 
a  single  time  counter  indicates  that  sufficient  collisions  have 
been  sampled  in  the  current  time  step  and  cell. 


6.4.1  Global  Collision  Time  Counter 

If  molecular  pairs  are  selected  for  collisions  such  that 
the  various  collision  classes  automatically  appear  with  the  proper 
relative  frequency  (see  below),  then  it  is  not  necessary  to  con¬ 
sider  separate  time  counters  for  all  the  various  collision  clas¬ 
ses.  One  approach  that  could  then  be  applied  is  to  keep  a 
collision  time  counter  for  just  one  collision  class,  and  incre¬ 
ment  it  when  collisions  of  that  class  occur.  If  the  various 
collision  classes  are  being  selected  according  to  their  correct 
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relative  frequency,  then  simulating  the  proper  frequency  for  one 
collision  class  will  ensure,  in  the  long  run,  that  all  collision 
classes  are  occurring  with  the  correct  frequency.  A  disadvantage 
with  this  approach  is  the  necessity  of  making  an  arbitrary  choice 
for  the  collision  class  which  is  to  have  a  time  counter.  Further¬ 
more,  there  may  be  no  good  choice  in  a  reacting  flow  where  the 
dominant  species  may  vary  strongly  from  place  to  place.  (Clearly, 
one  would  not  want  to  select  a  class  of  collision  that  does  not 
occur  in  a  given  cell,  since  the  result  would  be  a  never  ending 
sampling  of  collisions  of  other  classes.) 

The  preferred  approach  is  to  define  a  global  collision 
time  counter,  t^,  which  is  a  weighted  average  of  the  time  counters 
of  all  collision  classes;  i.e. 


i  E  “i:  "cij 

3=1 _ 


(72) 


where  the  are  nonnegative  coefficients  which  can  be  selected 
at  will.  Note  that  in  this  formulation  every  collision  will 
result  in  some  increment  of  the  global  time  counter  (unless 
=  0  for  that  class) ,  so  the  collision  sampling  frequency 
is  not  dependent  on  any  one  collision  class. 

It  remains,  of  course,  to  specify  the  ^ .  A  very  convenient 
choice  is  given  by 


n .  n . 

— i-i-r 
(1  +  6  .  . ) 
ID 


(73) 
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Firstly,  Eq.  (73)  is  convenient  because  it  tends  to  make  the 
collision  classes  with  the  higher  collision  frequencies  count 
more,  resulting  in  good  statistics  for  t  irrespective  of  cell 
location.  (Note  that  is  cell  dependent  since  the  species 
number  densities  <  e  cell  dependent.)  Secondly,  Eq.  (73) 
results  in  a  particularly  convenient  form  for  t^.  The  norm¬ 
alizing  factor  in  the  denominator  of  Eq.  (72)  can  be  summed  to 
give 


2 

T 

n 


E  E 

j=i 


n.n . 


(1  + 


Htt 


i=l 


13 


C13 


(74) 


Hence,  a  collision  of  class  i j ,  which  would  produce  an  increment 

of  At^. .  to  its  own  time  counter  produces  an  increment  At  to  t 
C13  g  g 

given  by 


At. 


CX3 


(75) 


where,  again,  n  is  the  total  number  density  of  all  species  in  the 
cell.  If  Eq.  (71)  is  substituted  into  Eq.  (75),  the  result  is 


At. 


Vn  a . . c 
13  r 


(76) 


Equation  (76)  is  extremely  significant  since  it  recaptures  the 
precise  form  of  the  time  counter  increment  for  a  single  species 
(Eq.  (64)),  but  indicates  that  it  is  completely  valid  for  a  multi- 
component  mixture  so  long  as  the  various  collision  classes  are 
sampled  with  the  proper  relative  frequency. 
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6.4.2  Collision  Pair  Selection  in  Multi-Component  Mixtures 

When  considering  selection  of  collision  pairs,  it  is  crucial 
to  remember  the  distinction  between  real  and  simulated  molecules 
discussed  in  Section  3.4.  Given  two  simulated  molecules  selected 
at  random  from  within  the  cell,  the  probability  of  their  having 
a  real  collision  is  proportional  to  W^W^a^jC^.  However,  real 
collisions  cannot  happen  individually;  they  come  at  a  time, 

ij 

where  W^.  is  the  lesser  of  W.  and  W..  Hence,  when  a  collision  is 

Li  1  3 

decided  upon  in  the  program,  W^  of  them  will  occur.  To  compen- 
sate  for  this,  potential  collision  pairs  should  be  accepted  for 
collision  according  to  the  size  of  Q  given  by 


Q 


"u“ij'=r 


(77) 


The  relative  frequency  of  real  ij  collisions  will  then  be  propor¬ 
tional  to  the  product  QW^^  (the  relative  probability  of  a  pair 
being  accepted  for  collision  times  the  number  of  real  collisions 
occurring  when  the  pair  is  accepted) ,  which  is  the  desired  rela¬ 
tion.  Selection  of  collision  pairs  with  the  correct  relative 
frequency  then  assures  that  incrementing  the  global  time  counter 
as  discussed  above  will  give  a  statistically  correct  sampling  of 
all  collision  classes  simultaneously. 


6.4.3  Summary  of  Collision  Sampling  in  Multi-Component  Mixtures 

The  results  of  this  section  can  be  summarized  via  the  fol¬ 
lowing  procedure  for  the  sampling  of  collisions: 

•  Each  cell  has  a  (current)  maximum  value  of  Q,  Q_  , 

max 

that  has  been  encountered  so  far  in  the  collision 
sampling  process.  Whenever  a  larger  value  is 
encountered,  is  set  equal  to  that  larger 

ItlaX 

value. 
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•  Each  cell  has  a  current  value  of  the  global  time 
counter,  t^. 

•  Pairs  are  selected  at  random  from  all  molecules 
within  the  cell. 

•  For  each  pair,  Q  (as  defined  by  Eq.  (77))  is 
computed . 

•  The  ratio  of  Q  to  is  computed,  and  a  random 

variable  is  generated.  The  pair  is  accepted  for 
collision  if  the  random  variable  is  less  than  the 
ratio.  (If  the  pair  is  not  accepted,  then  another 
random  pair  is  selected.  The  process  continues 
until  a  pair  is  accepted.) 

•  For  an  accepted  pair,  the  collision  mechanics  are 
computed  as  described  in  Section  4.  This  always 
corresponds  to  collisions. 

•  The  global  time  counter  is  incremented  by  W,.  At  , 

ii  g 

where  At  is  given  in  Eq.  (76). 
y 

•  The  process  continues  until  the  global  time 
counter  goes  beyond  the  overall  flow  time.  At 
that  point,  the  collision  sampling  is  commenced 
in  the  next  cell. 

•  When  all  cells  have  had  collisions  simulated,  then 
the  code  proceeds  to  the  translation  portion.  (See 
Sections  1  and  5.) 
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7.  INITIAL  AND  BOUNDARY  CONDITIONS 


7 . 1  General  Considerations 

The  initial  and  boundary  conditions  necessary  to  simulate 
a  problem  frequently  do  not  receive  their  fair  share  of  consid¬ 
eration.  It  is  these  conditions  which  usually  distinguish  one 
solution  from  another,  and  their  correct  and  efficient  specifica¬ 
tion  should  be  a  central  concern.  Nonetheless,  there  is  clearly 
room  for  advancements,  particularly  in  the  specification  of 
boundary  conditions.  Many  gas  dynamic  solutions  involve  boundary 
conditions  specified  at  infinity,  which  are  currently  simulated 
by  placing  boundaries  very  far  from  the  main  flow  region.  It 
would  result  in  a  substantial  computational  simplification  if 
boundary  conditions  applicable  closer  in  to  the  flow  region  of 
interest  could  be  generated  for  free  gas  boundaries.  Wall  bound¬ 
ary  conditions  also  frequently  involve  a  fair  degree  of  approxi¬ 
mation,  usually  taking  the  form  of  accommodation  coefficients 
(though  in  this  case  the  problem  is  as  much  a  lack  of  basic 
physical  understanding  of  the  gas-surface  interaction  process  as 
it  is  a  lack  of  a  good  ntimerical  simulation) . 


7.2  Initial  Conditions 

Since  the  direct  simulation  Monte  Carlo  method  is  inherently 
an  unsteady  technique,  an  initial  state  must  be  specified  in 
order  to  advance  the  solution.  (For  situations  where  a  steady 
state  result  is  desired,  it  is  obtained  as  the  long  time  solu¬ 
tion  to  an  unsteady  problem.  In  this  case  the  initial  conditions 
have  no  effect  on  the  eventual  solution,  but  they  may  well  have 
an  impact  on  the  speed  with  which  that  state  is  achieved.)  It 
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will  be  assiomed  here  that  the  initial  conditions  correspond  to 
a  uniform  flow  with  the  translational  and  internal  modes  being 
in  equilibrium.  The  specification  of  the  inital  conditions 
therefore  involves  determining  the  state  vector  elements  consis¬ 
tent  with  this  condition  for  the  desired  number  of  molecules. 

7.2.1  N\imber  of  Simulated  Molecules  and  Weighting  Factors 

The  desired  number  of  simulated  molecules  of  each  species 
in  each  cell  (referred  to  here  as  M  )  will  usually  be  an  input 
quantity.  (Typically,  simulations  aim  for  a  total  number  of 
molecules  per  cell  in  the  neighborhood  of  twenty.)  Given  the 
initial  number  density  to  be  simulated  for  a  species,  n^^,  (which 
will  have  been  converted  in  the  code  to  internal  dimensions  — 
see  Section  3)  the  weighting  factor  for  a  species  in  a  cell  can 
be  derived  from  an  application  of  Eq.  (26)  to  give 

=  Vn^/M^  ,  (78) 

where  V  is  the  cell  volume.  If  a  species  is  not  initially  pre¬ 
sent  in  a  cell,  then  the  weighting  factor  is  set  to  a  very 
small  but  positive  number,  typically  that  which  would  correspond 
to  an  initial  mole  fraction  of  one  part  per  million  or  so.  A 
weighting  factor  of  zero  would  cause  an  attempt  to  divide  by 
zero  when  molecules  of  that  species  move  into  the  cell,  but  it 
is  generally  best  to  start  the  weighting  factors  off  small.  As 
the  solution  proceeds,  the  weighting  factors  are  automatically 
adjusted,  but  the  adjustment  upward  is  more  direct  and  immediate 
than  the  adjustment  downward  (see  Sections  3.4  and  5.2). 
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7.2.2  Initial  Positions 

The  initial  molecules  assigned  to  a  cell  should  have  an 

equal  probability  of  being  placed  in  any  volime  element  of  the 

cell.  The  rules  for  accomplishing  this  will  change  with  the 

coordinate  system  being  used,  but  they  are  readily  derivable 

from  the  general  principle.  As  an  example,  consider  an  axisym- 

metric  simulation  where  a  cell  is  bounded  by  the  radial  positions 

r^^  and  r2  (r^  <  r2)/  and  the  axial  coordinates  2^  and  Z2  (z^  <  Z2). 

Since  the  basic  volume  element  in  this  coordinate  system  is 
2 

2TTrdrdz  =  irdCr  )d2,  the  volume  elements  will  be  sampled  with 

2 

equal  probability  if  r  and  z  are  sampled  randomly.  Hence,  a 
molecule  can  be  assigned  axial  and  radial  positions  via 

r  =  ^rj^^  +  3(r2^  -  (79) 

and 

2  =  z^  +  3(Z2  "■  z^)  .  (80) 

(Recall  that  every  time  the  symbol  g  occurs  it  implies  a  separate 
random  variable.  In  particular,  one  should  certainly  not  use 
the  same  random  variable  to  determine  radial  and  axial  coordin¬ 
ates.  The  second  application  would  hardly  qualify  as  "random".) 

7.2.3  Initial  Velocity  Components 

The  thermal  velocity  components  for  a  molecule  in  transla¬ 
tional  equilibrium  (neglecting,  for  the  moment,  any  mean  flow 
contribution)  should  be  selected  from  the  normalized  Maxwellian 
velocity  distribution,  f^ (v) ,  given  by 

fg  =  exp  [-(av)^]  ,  (81) 
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where 


“  "  yj  2R-T  ' 

m  is  the  species  molecular  weight,  Rq  is  the  universal  gas  con¬ 
stant  and  is  the  temperature.  Equation  (81)  applies  for  each 
of  the  molecular  velocity  components,  and  must  be  sampled  three 
times  for  each  molecule  that  comprises  the  initial  state  of  the 
simulation.  A  method  for  directly  sampling  from  this  distribu¬ 
tion,  as  presented  in  Ref.  1,  is 


^1 

=  2t:B  , 

(83) 

^2 

=  “^-log(B)/a  , 

(84) 

V 

=  A2sin{A^) 

(85) 

After  the  thermal  velocity  components  are  determined  for  each 
molecule,  then  any  mean  flow  velocity  is  simply  added  on.  The 
velocities  are  then  transformed  to  internal  units  (see  Section 
3.3) . 

7.2.4  Initial  Internal  Energies 

The  only  remaining  element  of  the  state  vector  to  be  speci¬ 
fied  is  the  internal  energy.  Internal  energies  for  a  gas  in 
equilibrium  are  distributed  according  to  the  normalized  distri¬ 
bution  function  f ^  given  by 

.C/2-1 

^i  ni72T 

where  ^  represents  the  number  of  internal  degrees  of  freedom  for 
the  species  in  question,  F  is  the  gamma  function  and  C  is  a 
dimensionless  internal  energy,  i.e. 


?  =  VV<»  ' 

where  is  the  internal  energy.  In  general,  sampling  of  Eg. 

(86)  must  be  done  via  the  acceptance-rejection  method.  In  the 
present  SSI  codes  z  is  restric-^ed  to  being  greater  than  or  equal 
to  two  (or  the  trivial  case  of  c  equal  to  zero,  which  just  gives 
Ej  identically  equal  to  zero).  -If  ?  is  precisely  equal  to  two, 
then  a  direct  sampling  of  the  internal  energy  is  possible  via 

C  =  -log(e)  .  (C  =  2)  (88) 

In  the  general  case  of  c;  >  2 ,  it  proves  convenient  to  first  intro¬ 
duce  the  transformation  s  =  VT".  s  is  then  distributed  in  propor¬ 
tion  to  the  distribution  g(s)  given  by 

g(s)  =  2s^^"^^  exp(-s^)  .  (89) 

Since  g(s)  is  tc  be  sampled  via  the  acceptance-rejection  method, 

xt  is  first  necessary  to  determine  its  maximum  value, 

Standard  calculus  serves  to  show  that  g  occurs  at  s  *  s*, 

max 

where 


s*  =  VU  -  l)/2 
so 


(90) 


g(s) 

g 

^max 


(91) 


The  sampling  of  Eq.  (89)  proceeds  as  follows: 

•  A  value  of  s  is  selected  randomly  from  the  interval 

s  .  to  s  ,  where 
min  max' 


max 


=  s*  +  5 


(92) 


and 


s  ■ 
min 


greater  of  {0,s*  -  5} 


(93) 
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(Note  that  in  this  interval  g(s)  goes  from  its 
maximum  value  to  a  value  on  the  order  of  10 
its  maximum  value.  The  transformation  from  Eq. 

(86)  to  Eq.  (89)  was  made  mainly  to  achieve  a 
probability  function  wnich  dies  off  extremely 
rapidly  away  from  its  maximum  value,  so  that 
very  little  error  is  associated  with  considering 
a  finite  subinterval  for  the  sampled  variable.) 

•  is  calculated  via  Eq.  (91). 

luaX 

•  A  random  variable  is  generated,  and  the  value  of 

s  is  kept  if  the  random  variable  is  less  than 

g(s)/g  .  Otherwise,  the  procedure  is  repeated 
in  3.x 

until  a  value  of  s  is  accepted. 

•  When  a  value  of  s  is  selected,  then  the  internal 

2 

energy  is  given  by  s 

As  for  all  the  initial  conditions,  the  codes  will  automatically 
then  express  the  values  in  internal  dimensions  (see  Section  3.3). 


7 . 3  Wall  Boundary  Conditions 

Boundary  conditions  in  the  direct  simulation  Monte  Carlo 
technique  are  applied  for  both  wall  and  free  gas  boundaries. 

In  a  wall  boundary  condition,  whenever  a  molecule  is  simulated 
to  strike  the  wall  some  action  must  be  taken.  Unless  the  wall 
is  a  condensing  boundary  the  molecule  will  be  reflected,  and  one 
of  several  boundary  conditions  can  be  selected.  The  easiest 
condition  is  that  of  a  specularly  reflecting  wall,  where  the 
molecule  simply  is  replaced  by  its  mirror  image  to  keep  it  within 
the  solution  region.  Another  frequently  used  condition  is  that 
of  a  perfectly  accommodating  wall.  In  this  condition,  a  molecule 


is  reemitted  from  the  wall  after  being  selected  from  a  distri¬ 
bution  characteristic  of  the  wall  temperature  and  velocity.  To 
date,  wall  boundary  conditions  have  not  been  required  in  the 
SSI  codes,  and  they  will  not  be  discussed  further  here.  They 
are  discussed  at  some  length  in  Ref.  1,  including  the  intermed¬ 
iate  case  of  partial  accommodation. 

7 . 4  Free  Gas  Boundary  Conditions 

In  many  cases,  the  boundary  condition  is  meant  to  simulate 
a  region  of  uniform  equilibrium  flow.  Molecules  leave  the 
solution  region  in  the  normal  course  of  their  trajectories, 
and  they  simply  disappear  from  the  simulation.  Molecules  are 
introduced  from  outside  the  boundary  as  selected  from  distri¬ 
butions  appropriate  to  incoming  molecules  in  the  undisturbed 
flow.  It  should  be  stressed  that  this  is  not  the  same  as  simply 
sampling  a  Maxwellian  velocity  distribution,  since  it  is  the 
molecular  flux  across  the  boundary  which  must  be  correctly  sim¬ 
ulated.  Hence,  molecules  with  a  large  component  of  velocity 
inward  from  the  boundary  are  more  likely  to  be  selected  than 
they  would  be  in  choosing  molecules  appropriate  to  a  static 
distribution  as  was  done  for  the  initial  conditions  described 
above . 

7.4.1  Incoming  Number  Flux 

The  first  requirement  is  to  determine  the  number  of  mole¬ 
cules  which  should  be  introduced  across  the  boundary  during 

the  solution  time  step  At  .  The  incoming  number  flux,  q, 

™  (1) 
(molecules  per  unit  area  per  unit  time)  can  be  expressed 

q  =  n^  I  exp  (-w^ ) /VT"  +  w  [l  +  erf  (w)]  |  /2a  ,  (94) 
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where 


w  =  au^cos(9)  ,  (95) 

and  6  represents  the  angle  between  the  inward  surfctce  normal  and 
the  mean  flow,  which  has  a  magnitude  (i.e.,  u^cos(6)  is  the 
inward  component  of  the  mean  flow  velocity),  n^  is  the  far  field 
number  density  of  the  species  of  interest  and  a  is  given  in 
Eq.  (82). 

Equation  (94)  must  be  applied  for  each  cell  on  the  boundary 
and  for  each  species  that  exists  in  the  ambient.  The  number  of 
simulated  molecules  to  be  introduced  into  the  cell,  Nj^,  is  given 
by 

Nj^  =  qA^Atyw  ,  (96) 

where  is  the  area  of  the  cell  along  the  boundary  and  W  is 
the  weighting  factor  for  the  species  and  cell  in  question. 

7.4.2  Incoming  Molecular  Velocity  Components 

A  coordinate  system  should  be  set  up  locally  at  the  bound¬ 
ary  such  that  one  direction  is  in  the  direction  of  the  inward 
normal  and  the  other  two  directions  are  perpendicular  to  it. 
Velocity  components  are  first  determined  in  terms  of  this  local 
coordinate  system,  and  then  transformed,  if  necessary,  to  the 
main  code  coordinate  system.  In  the  local  coordinate  system, 
the  velocity  components  parallel  to  the  surface  are  determined 
as  discussed  above  for  initial  state  molecules,  but  the  inward 
component  of  velocity  must  be  selected  in  proportion  to  the 
distribution  h(v)  given  by 

h(v)  =  (v  +  w)  [exp  -(av)^]  ,  (97) 
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which  must  be  sampled  via  the  acceptance-rejection  method. 
(Only  positive  v  is  allowed,  of  course,  from  the  definition 
of  the  coordinate  system. ) 


7.4.3  Incoming  Molecular  Position 

An  initial  entry  position  should  be  selected  for  the  mole¬ 
cule  such  that  the  flux  is  randomly  distributed.  Assuming  there 
is  no  variation  of  9  across  A  this  means  that  each  area  element 
of  the  exposed  cell  area  should  have  an  equal  probability  of 
being  selected.  Once  the  initial  entry  position  is  selected, 
then  the  molecule  should  be  translated  a  random  fraction  of 
At^  along  its  trajectory  to  determine  its  actual  location.  All 
of  the  considerations  discussed  in  Section  5  with  regard  to 
molecular  translations  apply  to  this  translation  and,  in  partic¬ 
ular,  it  must  be  possible  to  dynamically  adjust  the  weighting 
factor  as  required. 


7.4.4  Incoming  Molecular  Internal  Energies 

The  internal  energies  are  selected  as  described  in 
Section  7.2.4. 
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8.  STATISTICAL  SAMPLING  OF  OUTPUT 


8 . 1  General  Considerations 

It  is  safe  to  say  that  the  molecular  state  vectors  as  they 
exist  in  the  computer  do  not  comprise  the  desired  output  of 
the  procedure.  With  rare  exceptions,  it  is  usually  macroscopic 
quantities  such  as  temperature,  density,  mean  flow  velocity, 
etc.,  which  are  of  interest  —  not  the  microscopic  quantities 
represented  by  the  state  vector  components  of  an  individual 
simulated  molecule.  The  generation  of  the  desired  output 
requires  that  the  macroscopic  quantities  of  interest  be  repre¬ 
sented  in  terms  of  statistical  sums  of  the  available  microscopic 
quantities;  and  it  is  the  main  purpose  of  this  section  to  pre¬ 
sent  these  correspondences.  All  sums  are  kept  in  terms  of 
"real"  molceules  and  events,  i.e.,  the  current  weighting  fac¬ 
tors  are  included  in  the  sums.  This  is  essential  since  the 
weighting  factor  determines  the  statistical  importance  of  a 
given  molecule.  Since  the  weighting  factors  are  dynamically 
and  unpredictably  adjusted  as  the  solution  progresses,  it  would 
not  be  possible  to  go  back  and  add  in  the  effect  of  weighting 
factors  a  posteriori. 

In  general,  it  must  be  decided  ahead  of  time  exactly  what 
output  is  desired  from  the  code,  and  therefore  what  statistical 
sums  should  be  kept  to  generate  it.  There  is  a  vast  amount  of 
potential  information  in  the  simulation,  and  it  is  not  reasonable 
to  store  all  possibly  interesting  quantities  in  all  runs.  On 
the  other  hand,  it  is  wasteful  to  completely  rerun  a  case  just 
because  the  user  decides  there  was  an  additional  quantity  he 
was  interested  in.  The  selection  of  output  for  a  given  run, 
therefore,  unavoidably  requires  user  judgement.  Once  the 
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user  has  decided  upon  the  required  output,  the  determination  of 
which  statistical  sums  should  be  kept  can  be  done  automatically 
by  the  code.  Care  is  taken  to  make  sure  that  a  statistical  sum 
is  not  duplicated  internally  if  it  is  required  by  more  than  one 
requested  output  quantity. 

Some  initial  words  of  caution  are  required.  By  its  nature, 
the  direct  simulation  Monte  Carlo  method  works  with  far  fewer 
molecules  than  nature  does,  and  it  therefore  exhibits  consider- 
ably  greater  statistical  variation  in  its  macroscopic  predictions. 
To  reduce  these  variations,  the  codes  are  run  repeatedly  for  the 
same  case,  increasing  the  statistical  base  from  which  the  macro¬ 
scopic  output  is  derived.  If  care  is  taken  to  use  efficient 
techniques,  such  as  described  in  this  report,  then  useful 
results  can  usually  be  obtained  with  a  modest  computational 
effort.  This  statement  must  be  tempered,  however,  by  a  reali¬ 
zation  of  the  convergence  rate  for  Monte  Carlo  sampling.  Basic¬ 
ally,  the  statistical  error  in  the  output  converges  as  one  over 
the  square  root  of  the  sample  size  (or  run  time) .  Hence,  if  a 
solution  looks  good,  but  the  user  decides  he  would  like  one 
more  significant  digit  (i.e.,  he  would  like  the  statistical 
error  to  be  reduced  to  0.1  times  its  current  value)  it  would 
require  that  the  run  time  be  increased  by  a  factor  of  100!  It 
can  be  seen  that  the  desire  for  more  accuracy  can  very  quickly 
turn  the  most  efficient  code  into  a  money  gobbling  nightmare. 

When  using  a  Monte  Carlo  technique,  one  must  accept  some  sta¬ 
tistical  scatter  in  the  output. 

8.2  Sampling  of  Instantaneous  Output  Quantities 

Instantaneous  output  quantities  are  those  which  are,  in 
principle,  derivable  from  an  instantaneous  "snapshot”  of  the 
solution.  These  quantities,  such  as  density,  temperature  and 
velocity,  can  be  determined  by  examining  the  molecular  state 
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a 


vectors  at  a  particular  time  in  the  simulation.  The  code 
pauses  in  the  simulation  and  uses  the  molecular  state  vector 
elements  to  add  values  to  statistical  sxxms  appropriate  to  the 
various  cells  and  the  particular  time  that  it  paused.  It  then 
proceeds  with  the  simulation  until  the  next  sampling  time.  As 
the  code  goes  thrugh  its  successive  runs,  it  stops  at  the  same 
points  in  the  simulation  every  time  and  adds  to  the  statistical 
base  for  the  sums.  The  items  listed  below,  with  their  statis¬ 
tical  definitions,  are  selectable  as  output  requests  in  the  SSI 
codes.  Summations  are  made  over  all  applicable  simulated  mole¬ 
cules,  which  includes  separate  runs. 

Total  Number  Density 


Mean  Molecular  Weight 


(98) 


j'th  Velocity  Component 


(99) 


V. 

J 


W.m.V. . 
1  1  31 


W.m. 
1  1 


(100) 
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Overall  Translational  Temperature 


T 


where 


i 

’'i”!  fll  ''li  *  ''L) 

i 

®3  =  S 

i 

®4  = 

i 

S-  =  /  W.m.v_. 

5  1  1  3i 

i 

=6  =  H  "i-"! 


Translational  Temperature  in  j'th  Direction 


Tj  = 


^0^1 


^7  “  ^8/^6 


(101) 


(102) 


(103) 


(104) 


(105) 


(106) 


(107) 


(108) 
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where 


>_  =  ^  W.m.v?. 

7  1  1  Di 


So  =  >  W.m.v.. 

8  /  j  1  1 


Internal  Mode  Temperature 


22  ”i' 


■li 


'o  22  ”i'i 


(109) 


(110) 


(111) 


With  the  exception  of  Eq.  (99),  all  of  the  above  quantities  can 
also  be  defined  and  calculated  for  any  specified  species.  The 
sums  are  the  scune  except  that  only  molecules  of  that  species 
are  considered.  Before  printing  output  quantities,  they  are 
always  transformed  to  standard  dimensions  from  the  internal 
scales. 

8.3  Sampling  of  Time  Averaged  Quantities 

Some  additional  quantities  of  interest  are  not  sampled  at 
a  separate  sampling  time  as  described  above,  but  rather  as  the 
sir  lation  evolves.  Excunples  of  such  quantities  are  collision 
rates,  reaction  rates,  mean  velocities  between  molecules,  etc. 
In  general,  these  quantities  depend  on  the  relative  state  of 
more  than  one  type  of  molecule,  and  they  are  by  their  nature 


expressed  as  average  values  over  a  finite  time  interval.  The 
formulas  for  calculating  these  quantities  are  little  more  than 
event  counters,  and  will  not  be  included  here.  The  following 
quantities  are  currently  available  as  output: 

•  Mean  Relative  Velocity  Between  any  Two  Species 

•  R.M.S.  Deviation  of  Mean  Relative  Velocity 
Between  any  Two  Species 

•  Mean  Product  of  Cross  Section  Times  Relative 
Velocity  Between  any  Two  Species 

•  Collision  Rate  Between  any  Two  Species 

•  Reaction  Rate  for  any  Chemical  Reaction. 

The  sampling  for  all  of  these  quantities  occurs  in  the  collision 
simulation  routines.  As  pairs  are  considered  as  possible  colli¬ 
sion  partners,  statistics  are  kept,  if  necessary,  to  generate 
the  first  three  quantities.  Statistics  on  collisions  and 
reactions  are  kept  as  they  occur. 
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